%% hydraulic_calculation

load('data_PQT.mat');
P=zeros(size(d,1),length(lcgc));
T=P;
Z_s=P;
Z_s(:,1)=[Z Z Z]';
P(:,1)=[5 5 5]';
T(:,1)=[32 32 32]';
Q=[Q0 Q0 Q0 Q2 Q2 Q2];  %104Nm3/d
epsilonT=0.01;
epsilonP=10^-2;

for PlanNum = 1:size(d,1)
    P1=5;
    T1=32;
    Ppj=P1;
    Tpj=T1;
    for lc_x =1:length(lcgc)-1
        c=find(sectionPoint-lc_x>=0,1);
        d_x=d(PlanNum,c);
        D_x=D(PlanNum,c);
        thickness_x=thickness(PlanNum,c);
        Q_s=Q(c)/(24*3600);    %Nm3/s
        
        % first loop
        LoopSteps;
        T2_0=T2;
        P2_0=P2;
        for loopTimes=1:10
            LoopSteps;
            if (abs(T2_0-T2)<epsilonT && abs(P2_0-P2)<epsilonP)
                break;
            else
                T2_0=T2;
                P2_0=P2;
            end
        end     
        P(PlanNum,lc_x+1)=P2;
        T(PlanNum,lc_x+1)=T2;
        Z_s(PlanNum,lc_x+1)=Z_x;
        P1=P2;
        T1=T2;
        Ppj=P1;
        Tpj=T1;
        
    end
    
    
end
    Ppj_total=2/3*(P(:,1)+P(:,end).^2./(P(:,1)+P(:,end)));
    Tpj_tptal=[Tground(1) Tground(1) Tground(1)];